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Abstract 

We explore a computational approach to coarse graining the evolution of the large-scale features of a 
randomly forced Burgers equation in one spatial dimension. The long term evolution of the solution energy 
spectrum appears self-similar in time. We demonstrate coarse projective integration and coarse dynamic 
renormalization as tools that accelerate the extraction of macroscopic information (integration in time, self- 
similar shapes, and nontrivial dynamic exponents) from short bursts of appropriately initialized direct simu- 
lation. These procedures solve numerically an effective evolution equation for the energy spectrum without 
ever deriving this equation in closed form. 
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I. INTRODUCTION 



The behavior of physical systems is frequently observed, and modeled, at different levels of 
complexity. Laminar Newtonian fluid flow is a case in point: it is modeled at the atomistic level 
through molecular dynamics simulators; at a "mesoscopic" level through Lattice Boltzmann mod- 
els; and (in physical and engineering practice) through continuum macroscopic equations for den- 
sity, momentum and energy: the Navier-Stokes. What makes the latter, coarse grained description 
possible is an appropriate closure, in this case Newton's law of viscosity, which accurately models 
the effect of higher order, unmodeled quantities (such as the stresses) on the variables ("observ- 
ables") in terms of which the model is written (the velocity gradients). Such closures are often 
known from physical and engineering observation and practice long before their mathematical 
justification becomes available (in this case, through Chapman-Enskog expansions and kinetic 
theory). It is, of course, tempting to consider a situation in which the fine scale description is the 
Navier-Stokes equations themselves, and the coarse grained description is an evolution equation 
for some observables of turbulent velocity fields, or possibly of ensembles of such fields. Discov- 
ering the appropriate observables and deriving such evolution equations has been a long-standing 
ambition in both physics and mathematics. Our goal here is much more modest: we present a 
very simple illustration of a problem for which a fine scale description is available, and for which 
(based on extended observations of direct simulation) we have reason to believe that a coarse 
grained evolution equation exists - yet it is not available in closed form. We illustrate how (with a 
good guess of the appropriate coarse grained observables) we circumvent the derivation of a closed 
form coarse grained model, but still perform scientific computing tasks at the coarse grained level. 
The numerical tasks we demonstrate are coarse projective integration, which accelerate the (coarse 
grained) computations of the system evolution, and coarse dynamic renormalization, which (when 
the coarse grained evolution is self-similar, as the case appears to be here) targets the computa- 
tion of the self-similar solution shape and the corresponding exponents. This is an illustration of 
the so-called Equation-Free framework for coarse grained scientific computation in the absence of 



explicit quantitative closures and the resulting coarse grained evolution equations plL 
We consider the one dimensional in space, randomly forced Burgers equation, 

du Idu^ , 
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FIG. 1: Forcing amplitude 2^exp[~ ] and dissipation term Vj,ypfc 
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subject to periodic boundary conditions over the domain x E [0, In] . We start with random very low 
energy initial conditions. We are interested in the relatively long term and large scale properties 
of this system in the inviscid limit, when subject to a random forcing at the small scales. In 
order to achieve this inviscid limit numerically, we consider a general dissipation term of the 
form Vj^yp {—\)"^^d^"u/dx^" following Js]. The coefficients Vj^yp and n are chosen so that the 
dissipation acts only at the highest wavenumbers; here, their values are Vj^yp = 10^^^ and n = 1. 
This choice of dissipation is based on the assumption that the universal ranges of the system are 
not sensitive to a particular choice of the parameters Vj^yp and n [|30. The white-in-time random 
force f{x,t) is defined in Fourier space by 



m(o)fik',(D') = ^cxp 



{k-ksf 



5{k-k')5{(o-(o'), with^ 



(2) 



where k and (o are the spatial and temporal frequencies. The forcing term has a peak at wavenum- 
ber k = ks = 1000 and dies off at large wavenumbers; the dissipation term is essentially zero 
up to around wavenumber k = 3000 and only becomes important at much larger wavenumbers; 
see Fig. [B 

Before we start our computations, we briefly recall certain known features of the equations and 
their dynamics. It is quite remarkable that equations (|1I2|) have non-trivial asymptotic behavior 
in both the ultraviolet (UV, k > ks) and the infrared (IR, k < ks) limits. Attempts to use (|1I2|) 
as a simplified 1-D model for the investigation of small scale features of 3-D turbulence, aimed 
at recovering E{k) oc /:^5/3 (-j^g small scale limit, were based on deep similarities between the 
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two: the total energy in this system is ' /2 = K^{^Lfl^ and, as in 3-D turbulence, the inertial 
range <^k <^ where k^ is the dissipation wave number, is characterized by a constant energy 
flux Je = which, in the interval r <^L — l/ks, is reflected in the analytic behavior of the 
third-order structure function 



53(r) = (5,m)3 = {u{x + r)-u{x)y = -U^r. 



(3) 



Even in decaying high Re-turbulence this equation is correct: dE / dt is very small. At first glance, 
due to the shock formation leading to the ^(A;) °^ k^^ energy spectrum (rather than ^(fc) k^^ with 
the exponent x ~ 5/3), these attempts failed. However, the model revealed a non-trivial bi-scaling 



behavior of the structure functions S„{r) = {u{x + r) — r" for n < 1 and Sn{r) r for n > 

1 [I12I1 . The model also generated asymmetric probability densities with algebraically decaying 
tails P{dru) for large amplitude fluctuations drU < [13i|32|]. The properties of this tail attracted the 
attention of various groups and are reasonably well understood in terms of the dissipation anomaly 
introduced by [|25|] : they have also been analytically obtained in [12LI9|,110|]. The Burgers equation 



stirred by a force with the k~^ spectrum does lead to the logarithmically corrected Kolmogorov 
energy spectrum discovered in |4] and does not have multifractality Izst- 

The IR properties of (1 1121) in the range k <^ k^, on which we focus in this paper, are also 
quite interesting and non-trivial. It was determined, on the basis of the one-loop calculation (also 
numerically verified in 113311 ) that in the limit k <^ k^, (1 1121) is equivalent to the Burgers equation 



driven by a random force with f{k)-^ oc /c^-spectral density, which is called the KPZ equation, 
governing various interfacial phenomena jlTII . The renormalization group methodology applied 
to the KPZ equation lllin led to the equilibrium energy spectrum £(/:) oc ^ = const (here D = l) 
and to a non-Gaussian dynamic fixed point with non-trivial dynamic scaling exponents. 



coo^k^. 



(4) 



indicating strongly non-diffusive behavior: r{t) oc t2. These predictions, including dynamic scal- 
ing, have been numerically verified in ll33n . Numerical simulations in the IR range become 
progressively slower; our computer-assisted approach aims at accelerating such simulations (see 
also 1I20I1). 

Our numerical simulator integrates equation in time using a pseudo-spectral scheme with 
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a = 22000 Fourier mode spatial discretization. Thus, the smallest scale resolved by our simulator 
is Ajc = Tz/N = 1.428 x 10^^. The nonlinear term is computed using the fast Fourier transform. 
The temporal discretization treats the linear term exactly and uses a third order Runge-Kutta-Wray 



scheme for the nonlinear term. The random forcing term in Fourier space is f{k,t) = a/ ^/Ato^, 
where Ok is a random number generated from a Gaussian distribution with zero mean and unit 
standard deviation, and At is the time step of integration. The various parameters chosen are At = 
5 X 10"^, /o = 5 X 10^, ks = 1000, Gf = 250, Vj^yp = 10"^^^ ^nd n = 7. An ensemble of 16 distinct 
realizations was computed, and the results averaged to achieve good statistics. 

We begin by observing direct numerical simulations of the evolution of equation ^ over time. 
Fig. [2] shows such a transient, starting from the zero initial condition, u{k,t = 0) = 0, for all k. The 
energy spectra ^(fc,?) are plotted after 1000, 5000, 25000, and 300000 time steps; the insets show 
the corresponding ID velocity fields. Visually, we can discern no clear trend in the evolution of 
the velocity fields themselves, which appear random; yet, we easily see a definite smooth structure 
in the evolution of the energy spectrum. 

For k > 700, the spectrum seems to quickly achieve a steady state profile within less than 
1000 time steps: the energy provided by the noise is balanced by the strong dissipation at high 
wavenumbers. It is much more interesting to observe E{k) for I < k < 500; here the spectrum 
quickly acquires a "corner-like" shape which appears to evolve smoothly in time, and moves to- 
wards lower wavenumbers while maintaining this shape. This "traveling front-like" evolution 
progressively slows down as the "comer" evolves towards lower wavenumbers. 



Theoretical considerations III , 



33|] indicate that in the range k <ks, the evolving energy spec- 



trum E{k) consists of three regimes. For k = 0{ks), the spectrum is not universal and depends on 
the properties of the forcing function. For l/L<^kj{t) <^k<^ k^, E{k) = k^^^ = const (here D = 
1); and in the time-dependent large-scale limit 1/L <^ k < kj{t), the energy spectrum is given by 
a universal relation E{k) oc k^~^^^ oc which we observe in our simulations. Superficially, the 
temporal evolution we observe resembles the time-dependence of a two dimensional flow driven 
by a small {ks = 0{\)) random force, eventually leading to the energy spectrum piling up on the 
smallest (integral) wave number ki = n/L and leading to the IR asymptotics E{k) oc k^^. We 
stress that the situation we are interested in may be quite different: in 2D, due to enstrophy con- 
servation in the infrared range (k < ks), there exists a constant, nonzero, energy flux toward large 
scales, leading to the so called "inverse cascade". The third-order structure function in a 2D flow 
is 53 = i^r 7^ 0. In the one dimensional flow we are interested in here, the large scale ther- 
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FIG. 2: Evolution of E{k) from the zero initial condition u{k,0) = 0, for all k. The spectra are shown 
after 1000, 5000, 25000, and 300000 time steps, with time increasing in the direction of the aiTow. The 
velocity fields coiTcsponding to the spectra are shown in the insets, in the same color as the spectra. 



modynamic equilibrium is flux-free and ^3 must be equal to zero, provided the system is large 
enough, that is, kj /ks <ti 1. 

We do not work at this very long time regime; we rather observe the evolution at times inter- 
mediate between this regime and the (relatively) short time it takes for the spectrum to stabilize 
for k > 700. This is the regime in which smooth and, as we argue below, apparently self-similar 
evolution of the energy spectrum prevails. 

Based on the study of many transients like those in Fig. [21 we decided to observe the system 
dynamics in terms of the energy spectrum E{k,t) only, which is defined as 

E{k,t) = {u{k,t)u*{k,t)), (5) 

where u{kj) is the spatial Fourier transform of the velocity field u{xj), and (■) denotes an ensem- 
ble average (here over the 16 realizations). A formal diagrammatic expansion for the evolution 



of the energy spectrum corresponding to (|1I2I) can be developed yOf l . In the limit k ^ 0, taking 
into account that all odd-order correlation functions rapidly tend to zero, the solution can be ac- 
curately represented in terms of an infinite series involving complicated integral expressions with 
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(a) 



(b) 





FIG. 3: Long term solution, after 1.25 x 10^ time steps (a) E{k) vs. k; (b) Si{r) vs. r/Ax, where Ax = n/N 
is the smallest scale resolved in our simulations. 



kernels which are convolutions of the energy spectrum E{k) only. Based on this, we assume that 
an evolution equation of the general form 



dE{k,t) 
dt 



(6) 



exists and closes for the (ensemble averaged) E{k,t). We do not know the nature of the operator ^ 
- it may not be local (in fc-space), and the equation may not be a partial differential equation - 
yet the implication is that knowing E{k,0) is enough to predict the (realization ensemble aver- 
aged) E{kj) for later times. If an evolution equation closes with respect to the expected spectrum 
only, then the remaining degrees of freedom (the phases of the velocity field Fourier coefficients) 
do not explicitly appear in this equation. This implies that either these variables do not couple to 
the spectrum evolution, or that their expected dynamics become quickly slaved to the spectrum 
evolution (with a relaxation time that is small compared to the spectrum evolution time scales). 
We test this latter ansatz by observing the solution ensemble not only through the spectrum t) 
but also through the evolving third order structure function 



S3(r,0 = Uu{x,t)-u{x + r,t)f 



(7) 



For completeness, the "extremely long term" steady state of E{k,t) and S3(r,?), obtained after 
integrating for 1.25 x 10^ time steps is included in Fig. [3l Fig. [3ta) suggests that E{k) reaches 
a constant value of 2 x 10^"^ in the range 1 < ^ < 150. The structure function ST,{r), plotted in 



7 



Fig. [2b), is practically zero in the range 200Ax; <r< 2k. The steady state value of S^{r) is nonzero 
for r < 200 Ax, due to the force acting primarily at small scales (corresponding to k = 1000) and 
strong dissipation at the smallest scales. 

The remainder of the paper is organized as follows: In Section HI] we discuss simulations sug- 
gesting an apparent dynamic self-similarity in the intermediate-time evolution of the spectrum. 
Section Hn] briefly outlines equation free computational techniques, and demonstrates coarse pro- 
jective integration of the unavailable equation for the expected spectrum evolution. Section |IV] 
then shows how to implement a dynamic renormalization fixed point algorithm that converges on 
the self-similar shape and pinpoints the temporal similarity exponent of its evolution. We conclude 
with a brief discussion. 

II. APPARENT DYNAMIC SELF-SIMILARITY 

Here we argue that the transient evolution of E{k) exhibits dynamic self-similarity. Fig. |4] 
clearly shows that the qualitative shape of the spectrum in the region I <k < 500 remains effec- 
tively unchanged, while gradually "travelling" to the left, slowly filling up the lower wavenumbers. 
We make this observation more precise as follows: We define a "fulcrum point" k = kf as the left- 
most end of the apparently steady region of E{k); here kf ^ 450. As shown in Fig. Ufa), we 
draw straight line segments of various slopes my (logarithmic in E as well as in k) starting at the 
fulcrum and intersecting successive computed spectra. Let dmj (t) be the distance of the point of 
intersection of each such line with the spectrum at time t. Our ansatz is that the distance d,nj{t) 
becomes uniformly stretched in (logarithmic) time. This is tested as follows: let the stretching 
factor along a line of slope mj be rmjit) = d,nj{t)/d,nj{ti), where ti is a reference time. Fig.Hfb) 
is a plot of r,„j versus my at various times in the order listed in the figure caption; it strongly sug- 
gests that the spectrum gets stretched by the same amount along its entire profile (notice that no 
intersections arise with my < —0. 1). Fig.Hfc) shows the time evolution of the stretching factor r{t) 
(averaged over the lines of different slopes), clearly showing linearity in logarithmic time. Thus, 
the transient evolution of E{k) appears dynamically self-similar: the spectrum in the wavenumber 
region 1 < ^ < 450 becomes uniformly stretched in logarithmic time, while maintaining its shape. 
We take advantage of this information in speeding up computations in the next section. 
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(a) (b) (c) 




\ogk slope, m \og{t/At) 



FIG. 4: (a) Schematic diagram showing calculation of the stretching factor. For a hne of slope m, r,„{tn) = 
dm{tn) /d,n{t\)- (b) Stretching factor r„,(?,) vs. slope m after 7, 10,25,50, 100, and 300(x 1000) time steps; 
the direction of time is indicated by the arrow. Stretching is measured with respect to the spectrum at 
t = 2500A?. (c) Mean stretch r{t) vs. log(?/A?) (solid line) and a least squares fit (dotted line). 



III. COARSE PROJECTIVE INTEGRATION 

Based on our simulation observations above, we now hypothesize that an evolution equation 
exists and closes at the level of the energy spectrum E{k,t) averaged over our ensemble of real- 
izations. If this equation were explicitly available, a typical numerical integration scheme would 
start by discretizing it in k (for example through finite differences, FD, or finite elements, FE) and 
then proceed to the derivation of ordinary differential equations (ODEs) in time for the vector of 
values a at the FD discretization points, or of the FE coefficients. These ODEs would be of the 
general form dsi/dt = F(a); the temporal evolution would be obtained using established numerical 
integration techniques, of which the simplest is arguably the explicit forward Euler method 



a(?o + AO =a(?o)+A?F(a(?o)). 



(8) 



Clearly, implementing this, as well as other initial value solvers, requires an evaluation of the time 
derivative (the right hand side of the ODE set) at prescribed system states. In projective inte- 



gration [[13 



1411 this time derivative is not obtained through such a function evaluation of F(a), 



since F(a) is not available in closed form; instead, it is obtained from processing the results 
of short bursts of simulation: starting at time tQ, we integrate the original equation for a short 
time with a time step 5t <^ A? to obtain a sequence a(ro),a(ro + 5?),a(?o + 25?), . . . , a.{tQ-\-n5t), 
where ndt < At. We can then use, for example, the last two points of this sequence to estimate the 



9 



time derivative, and predict a{to + At) as 

a{to + At) ^ a{to + {n-l)5t) + — — '^^^ ^''^^ {a{tQ + n5t) - a{tQ + (n - 1) 5?)) ; (9) 

better estimators (least squares, maximum likelihood) of the time deerivative can clearly be used as 
well. The procedure provides an approximation of the "projected in time" state a(fo + Af), which 
can then be used to restart the computation, and the process repeats. 

When the only available tool is a "black box" simulator of the system evolution equation, 
projective integration may significantly accelerate computations for systems with gaps in their 
eigenvalue spectra (separation of time scales). However, in our case, the equation we want to 
projectively integrate is an equation for energy spectrum evolution, while the available simulator 
is a simulator for velocity fields; we want to projectively integrate the coarse grained behavior 
(spectra) through observations of the fine scale behavior (velocities). This now becomes coarse 
projective integration. From our simulation we have available the full state (the velocity field 
ensembles) at time to; full direct simulation will then provide full states at each of the intermediate 
time steps. We restrict the full state to the coarse observables (the discretization of the spectrum); 
we use these observations to estimate the time derivative of these coarse observables; and we 
finally pass the time derivative estimates to our integrator of choice (here a simple forward Euler) 
in order to produce our approximation of the coarse observables at the next (projected) point in 
time to + At. 

Here comes a vital step in our equation-free coarse graining approach: our "black box" simula- 
tor requires a fine scale initialization (detailed velocities); yet we only have coarse grained initial 
conditions (ensemble-averaged spectrum discretizations). For the procedure to continue, we need 
to somehow construct full fine- scale states consistent with the desired coarse initial conditions. 
This construction is the so-called lifting step [|16L The missing degrees of freedom in 



our construction are the phases of the velocity fields. Our assumption that an equation exists and 
closes in terms of only ensemble averaged spectra implies that these remaining degrees of freedom 
either do not couple with the evolution of the coarse grained variables we retain, or, if they do (as 
we expect), they quickly get slaved to the coarse variables. "Quickly" here implies a compari- 
son with the time scales of the evolution of the retained coarse variables themselves. Exploiting 
this "quick slaving" (which is implicit in the assumption that an equation for the retained coarse 
variables exists and closes) is the main point of equation-free computation: since the fine scale 
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code is available, running it for a short time, even with wrongly initialized "additional variables", 
should provide the required slaving on demand. This can be considered a subcase of optimal pre- 
diction IlZI]. Given an energy spectrum initial condition, we therefore initialize the velocity fields as 
follows: the amplitudes of the velocity Fourier modes are chosen consistently with the prescribed 
initial energy spectrum, while the phases are chosen randomly. We then let our direct simulator 
evolve for some initial "healing" interval (described below) before we start observing the solution 
energy spectrum, in order to estimate its time derivative; this initial period is allowed so that our 
"wrong" phase initialization becomes "healed", that is, so that our phases or their statistics be- 
come slaved to the evolution of the spectra. We monitor this slaving through observations of the 
structure function S^irj) defined earlier. 

Here is a concise description of our coarse (spectrum level) projective (forward Euler) integra- 
tion. The procedure starts after we have already simulated (starting with zero initial conditions) 
for a period of approximately 3000 time steps; during this period the "high-fc" end of the spec- 
trum ^(fc > kf) approaches stationarity. 

1. Starting with a "current" velocity field initial condition, we evolve the randomly forced 
Burgers simulator for a short interval T, and restrict the velocity fields at successive in- 
stances within this interval to ensemble-averaged energy spectra. Fig. |2a) shows the re- 
stricted spectra at times t/At = 3500, 4500, 5500, and 7000. 

2. Using the procedure described in section [III we compute the stretching factors r,„j{t) along 
various straight lines passing through the fulcrum and intersecting these spectra. Our dis- 
cretization of the spectrum does not come in the form E{ki) for a number of mesh points hi 
(or log/:,), but rather in the form {E{mj) , k{mj)) for a number of "discretization slopes" my, 
the ensemble-averaged spectrum evolution appears smooth enough, so that approximately 
20 such discretization angles and simple interpolation is sufficient. We now compute the 
rate of change of the stretching factors in logarithmic time, and use the same to project the 
energy spectrum to a future time ^ T . Here, the stretching factors r„j . are computed along 
straight lines of slopes —0.2 < < 0.1, and their rate of change in logarithmic time is 
computed using a least squares fit as 

r^.{t)^P";^+P'^'\ogt, (10) 
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where ^ and ' the coefficients of the least squares fit. We then use equation (UOl) to 
project the stretching factor to a time t* ^ T, and compute the projected spectrum using 

\og[E{mj, f)] = r,„.(?*)log[£(mj, 
and \og[k{mj,t*)]=ryn.{t*)\o%[k{mj,ti)]. (11) 

An estimate of the spectrum at t* = SOOOOA? using those at times f/A? = 3500, 4500, 5500, 
and 7000 is shown in Fig. Oa). 

3. We now lift the projected spectrum to an ensemble of consistent velocity fields in order 
to repeat the process. This is done by randomizing the phases for all the wavenumbers, 
while retaining the amplitudes corresponding to the projected spectrum. More precisely, 
we generate 16 different initial velocity fields as u{k,t) = ^/E{kj)exp{i^k), where is a 
random number between and 2n chosen from a uniform distribution. 

4. We now repeat the process: we evolve the fine scale simulator (the randomly forced Burgers) 
with the initial conditions thus obtained; but before we start using the resulting spectra to 
estimate a new coarse time derivative and make a new forward Euler projection, we test the 
slaving of the injected degrees of freedom (the random phases) to our coarse observables 
(the spectrum). This slaving is monitored through the evolution of the newly initialized 
structure function 53(r). Fig. [3b) shows the structure function from a simulation that has 
been running for some time; if we restrict to energy spectra, and then randomize the phases 
(that is, lift back to velocity fields) this phase information is lost. When we restart the 
simulation from the lifted velocity fields, however, we clearly observe that the phases "heal", 
that is, 53 acquires its original shape quickly compared to the evolution time scales of the 
spectrum (in less than 500 time steps, which is much smaller than the projection horizon of 
> 20000 time steps). We continue now the "short burst" of fine scale simulation to collect 
data (restrictions to spectra) for the next projection step; the procedure (steps 1-3) then 
repeats. 

Fig. |5] shows an instantiation of such a coarse projective step. The restriction of an initial simula- 
tion used to estimate the coarse time derivative are obtained att/At = 3500, 4500, 5500, and 7000; 
the projection step (which here is performed in logarithmic time) brings us to ? = 30000A?; the 
projected spectrum (black dashed line) is compared to the result of continuing the full simulation 
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(a) (b) 




1 2 3 4 5 6 7 10° 10' 10' 10' 10' 



log A: r/Ax 

FIG. 5: (a) Coarse projective integration and its comparison with the original direct simulation (spectra from 
the continuing original simulation are shown in blue). The phases of all the modes are randomized after 
projection, (b) Healing of ST,{r) (black) in 400 time steps. After randomizing the phases att = SOOOOAf (flat 
line), their transient recovery is depicted after time steps 10, 20, 30, 40, 50, and 400. Also shown is a 
comparison with the original simulation att = 30400 (blue). 



until time t = 30000A? (blue). More importantly, the result of lifting from the projected spectrum 
and continuing the evolution till t = 70000A? (black dashed line) is also seen to coincide with 
the result of a full simulation till t = VOOOOAf (blue). Also, the computational savings in the first 
step of coarse projective integration can now be quantified. Using coarse projective integration, 
the randomly forced Burgers equation is simulated for 7000 time steps and then a projection step 
is used to obtain the spectrum at t = 30000A?, thus resulting in a saving of 23000 time steps, or 
a speedup factor of 4.3. The factor is actually a little lower because of the computational effort 
required for the projection step. Since the dynamics evolve in logarithmic time, this factor will 
increase for the later projective integration steps. 

The choice of performing projective integration (in effect, Taylor series expansion of the coarse 
solution in the independent variable) in logarithmic rather than linear time is motivated by the ob- 
servation that our transient is close to a dynamically self-similar solution. A detailed discussion of 
the usefulness of of (coarse) projective integration in problems with continuous symmetries, such 
as scale invariance, having (possible approximately) self-similar solutions can be found in jlSl] . 
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IV. COARSE DYNAMIC RENORMALIZATION 



Direct simulations starting at E{k) =0 showed that the spectrum quickly evolves to a "steady 
state" for k > kf ^ 450; after this initial period, its low wavenumber portion acquires a "constant 
shape" which is then observed, over time, to stretch gradually and rather smoothly towards k = 0. 
This stretching becomes increasingly slower in linear time. In what follows we do not concern 
ourselves with the very fast initial transient (the establishment of a steady state for k > kf); nor 
will we study the very long term dynamics that occur when the "corner" of the spectrum finally 
approaches ^ = 0. We focus on the "intermediate time" regime during which the k < kf spectrum 
acquires its "constant shape" which then appears to travel (in logarithmic k and logarithmic time) 
towards k = 0. This latter observation suggests an even simpler caricature of the self-similar 
evolution; see Fig. 13 a): the spectrum to the left of kf acquires quickly the appropriate "comer 
like" shape, which we approximate by two straight lines meeting at a point with wavenumber (0 • 
The energy of the wavenumbers to the right of ko remains more or less unchanged, anchored at the 
fulcrum point. The portion of the spectrum to the left of ko{t) subsequently simply moves towards 
lower k with constant shape. A constant shift of this part of the spectrum in log/c corresponds 
to a stretching by a constant factor in k, and clearly suggests that for our apparently self-similar 
solution 

E{k,t)=E(^k{t + tof^, for k<ko{t), (12) 

and an appropriate choice of constants /3 and to. For the data shown in Fig. |4l initialized with a 
zero spectrum at f = and taken at times t > 125A?, the values of the exponent and the constant to 
are computationally found to be /3 ~ 0.633 and to ~ —9SAt. These values can be extracted from 
Fig-Etb), which shows plots of the wavenumbers ko{t) versus t/At and ko{t) versus {t + to)/At on 
a logarithmic scale. 

This appears consistent with the numerical simulations of Yakhot and She [jssl as well as with 
dynamic renormalization group results suggesting a non-trivial dynamic scaling of the temporal 
frequency (O = 0{k^/'^) (or kt^/^ ^ const) as opposed io (O °^ vk^ or k = O((0^/^). The same 
exponent is obtained using different spectrum "caricatures" based on the same low-wavenumber 
straight-line approximations (not shown here). We would like to stress that computation of dy- 
namic scaling exponents is an extremely difficult and computationally expensive task. The fact 
that the non-trivial exponent Q emerged from our relatively fast calculation seems remarkable. 
Moreover, it serves as an independent test of our coarse graining procedure. 
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(a) (b) 




FIG. 6: (a) A cartoon illustrating the approximations to the spectra used to compute the self-similar solution. 
The plot shows spectra at a few select time instances. Their straight line approximations used to compute 
the exponents are shown in red. (b) Wavenumbers at the intersections of the spectra with the dotted green 
line are used to fit the constants j8 and to (see text). A plot of these wavenumbers vs. t/At is shown (red, 
o), with a least squares fit of the form ko{t) = ^const(^ + ^0)^ (black, □). Also shown are plots of ko 
vs. {t + to)/At (blue, V) and the corresponding least squares fit (black, +). 



Given these observations, it is clear that the intermediate time behavior can be recovered in an 
even more economical fashion than coarse projective integration: it is enough to find the "right 
shape" (the slope of the leftmost part of the spectrum) as well as the correct stretching rate (en- 
coded in the exponent); this information suffices to approximate the spectra at all intermediate 
times. This is analogous to computing the shape and the (constant) speed of a traveling wave 
in a problem with translational invariance. Knowing the shape and speed we can construct the 
solution at any future time; scale invariance here takes the place of translational invariance for 
traveling, and the stretching rate (quantified by the similarity exponents) is analogous to the trav- 
eling speed. The problem therefore reduces to finding the "right" spectrum shape (for k < kf): a 
spectrum E{k,t\) which, when evolved for some time T to E{k,t\ + T) and rescaled in /c-space, 
remains unchanged. This gives rise to the fixed point problem: 

£ft<,)-£(*(^ii^)'.n+r)=o. (13) 

Solving such coarse fixed point problems using coarse time steppers has been discussed and il- 



lustrated in several contexts 



ointproDJ 



3411 ■ The simplest fixed point algorithm is successive sub- 



stitution: we start with an initial guess of the self-similar spectrum shape which we evolve for 
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some time using the direct randomly forced Burgers simulator, including ensemble averaging. We 
then rescale the resulting spectrum by a constant factor, and repeat the process; for a discussion 



of a so-ca^ 
see 



led template based approach to choosing the appropriate scaling factor at each iteration, 
270, and our choice for this problem will be discussed below. If the self-similar spectrum 
is stable (as is the case here), this iteration converges to its scale-invariant shape. We remark, 
before proceeding, that other fixed point algorithms (in particular, matrix-free Newton-Krylov- 
GMRES, H) have been used for solving this type of problem; such Newton-based algorithms 
can converge to self-similar solutions even if they are not stable. In what follows, we demon- 
strate that successive substitution using a dynamically rescaled coarse time stepper does indeed 
converge to the right (invariant) spectrum shape. 

Long direct simulation will also eventually converge to this invariant shape; yet, while ap- 
proaching the right shape, the transient will also move towards lower wavenumbers (will "stretch") 
and this makes the evolution increasingly slower. We therefore have the potential to realize a com- 
putational advantage by rescaling the spectrum after partial evolution to a constant reference scale, 
a scale in whose neighborhood the dynamics are as fast as possible in physical time; this allows 
us to observe the approach to the stable invariant shape without this shape stretching. 

We now consider another caricature of the (ensemble averaged) spectra during our 
intermediate-time simulation, approximated by two straight lines in the region k < kf as, shown in 
Fig. Ha). These two straight lines intersect at a "corner"; we call this comer wavenumber k = kc. 
As we discussed earlier in this section, we observe that the slope of the straight line approximating 
the region of the spectrum to the left of the comer remains (approximately) constant (m = m^^f^ 
during self-similar evolution; finding the right shape becomes then equivalent to finding this slope. 
In order to obtain the self-similar shape, we initialize our spectrum as follows: The portion of 
the spectmm to the right of the fulcrum E{k > kf) is initialized at its (already) steady shape. 
For kc < k < kf, a straight line approximation is used, the slope of which is computed by a local 
fitting to the spectrum obtained from the initial direct simulation. Finally, for the region k < kc, 
we initialize the spectrum with a straight line of a "guess" slope m^Qf^. We now pin our refer- 
ence scale by selecting a particular, fixed kc-; this means that our procedure will converge on the 
representative of the one-parameter family of self-similar shapes that has its comer fixed there. 

Our coarse renormalized time stepper involves lifting from this guess spectrum to full Burgers 
velocity fields, (ensemble) evolution of these fields by the direct randomly forced Burgers simu- 
lator for a time T, restriction back to spectra, and ensemble averaging, after this time. A further 
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restriction step is now available by our choice of caricature, in the form of approximating the spec- 
trum by two piecewise linear segments. This gives us a new corner kc and a new slope tn^Qf^. We 
now rescale the resulting spectrum by shifting back to the old kc but keeping the new m^^j-f, this 
allows us to persistently observe the system at a convenient scale. Iterations of this "lift-evolve- 
restrict-rescale" procedure is shown in Fig. Ha); Fig.|7];b) shows a sequence of several successive 
renormalization substitutions converging to the invariant shape; here it takes approximately six 
iterations to convergence if the originally guessed slope is less steep than the true one; steeper 
initial slopes practically converge within one iteration. 

From the last iteration (upon convergence) we obtain the shape of the leftmost part of the spec- 
trum E{k,t) ^i-93±o.05 m^Qf^ ^ 2; see [jl IQ). We can also extract the scaling exponent by 
doing a similar least squares fit to a plot of log^o versus log?. The constants in the expression (fT2l) 
are /3 = 0.635 and tQ ^ —6685 A? (here we started iterating dXt= lOOOOAf of the original simu- 
lation). The larger we can make the comer kc (the closer to kf we can practically choose it), the 
fastest the approach to the self-similar shape (the dynamics slow down significantly at progres- 
sively lower wavenumbers). 
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V. SUMMARY AND DISCUSSION 



Using a particular type of randomly forced Burgers equation in one spatial dimension as our 
illustrative example, and based on direct simulations, we postulated that the system dynamics 
could, under certain conditions, be coarse grained to an effective evolution equation for the system 
energy spectrum. The direct simulations also strongly suggest an effectively self- similar evolution 
of this energy spectrum; it is important to reiterate that this ansatz is only made for a particular 
(intermediate) wavenumber regime, and only for initial conditions that contain no appreciable 
energy in the large scales. Based on these observations we demonstrated the implementation and 
use of the "equation-free" computational methodology: the design of short bursts of appropriately 
initialized computational experiments with the direct simulator that accelerate the study of the 
system behavior - as compared to full direct simulation only. In particular, we demonstrated 
coarse projective integration, as well as coarse dynamic renormalization. 

Our computations are predicated upon advance knowledge of the right set of coarse observ- 
ables; these are the quantities in terms of which the (unavailable) evolution equation can in prin- 
ciple be written. Our coarse observables here consisted of (a discretization of) the (ensemble 
averaged) energy spectrum; the phase variables were assumed (and observed, through the use of 
the third-order structure function 53) to become relatively quickly slaved to the energy spectrum 
dynamics. Algorithms that attempt to detect such observables from data-mining simulation results, 
such as the diffusion map approach Isl, 24], may prove valuable in extending the applicability of 
the computational tools we presented here. The computational construction of synthetic turbulence 
fields consistent with such observables (see, for example, the recent paper 112611 ) can be considered 
as a type of "lifting" scheme in our framework. Beyond the design of (numerical) experiments 
to accelerate system modeling, it is also possible to use such experiments to answer questions 
about the nature of the unavailable equation (for example, whether it is a conservation law). Our 
approach made no assumption about the local or nonlocal nature of the (explicitly unavailable) 
evolution equation: it may be a partial differential equation or a nonlocal, integro-differential one. 
It would be particularly interesting to devise computational experiments with the direct simulation 
code to explore whether the right-hand- side of our unavailable evolution equation appears to be 
local or not in A;-wavenumber space. 

The method proposed here may prove useful for the efficient extraction of dynamic scaling 
exponents in models of various physical phenomena for which we do not know the governing 
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equations in closed form. 
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